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CHAPTER  1 

INTRODUCTION 

Multistage  optimization  problems  of  management  systems  arise  in 
connection  with  processes  developing  in  time  in  which  one  or  more  control 
variables  must  be  controlled  to  achieve  certain  conditions.   Out  of  all 
possible  values  for  the  control  variables,  the  one  which  Rives  a  certain 
maximum  (or  minimum)  performance  index  while  simultaneously  keeping  all 
the  state  and  control  variables  within  the  specified  constraints  of  the 
problem  must  be  determined. 

Various  methods  have  been  developed  and  applied  to  solve  multistage 
and  single  stage  problems.   The  Kuhn-Tucker  method,  quadratic  programing 
and  linear  programming  is  in'  a  sense  an  adjacent  extreme  point  method 
which  employs  the  simplex  algorithm  as  the  fundamental  computational 

tool. 

The  calculus  of  variation  is  a  classical  tool  for  solving  optimization 
problems.  Until  recently,  however  computational  difficulties  limited  this  to 
solving  simple  problems  only.  Dynamic  programming  provides  an  entirely 
new  concept  of  optimization  and  has  been  used  quite  extensively  to  solve 
management  optimization  problems.  However,  owing  to  the  dimensionality 
difficulty  and  the  limited  fast  memory  capacity  of  present  day  computers , 
this  technique  cannot  be  applied  to  problems  with  a  fairly  large  number 
of  state  variables.   In  view  of  the  complexity  of  many  industrial  and 
management  systems,  this  is  a  serious  limitation. 

The  gradient  technique  or  the  method  of  steepest  ascent  is  an 
elementary  concept  in  solving  optimization  problems.  It  dates  back  to 


Cauchey  [Ita]  and,  in  a  variational  version,  to  Hardomarrl  [5a].   It  is 
based  on  the  fact  that  if  movement  is  made  in  the  direction  of  the  gradient 
of  the  objective  function,  movement  is  also  made  in  the  direction  of 
the  maximum  rate  of  increase  in  the  objective  function.  In  recent  years 
the  computational  appeal  of  the  method  has  led  to  its  adaptation  in  a 
variety  of  applications.  The  dynamic  version  of  the  technique,  generally 
known  as  the  functional  or  serial  gradient  technique,  has  been  applied 
successfully  to  solve  problems  in  aerospace,  control  and  chemical 
engineering  systems,  [lt-6,  8-13] 

This  report  is  a  study  of  the  way  in  which  the  gradient  concept  can 
be  applied  to  the  solution  of  optimization  problems  with  constraints. 
The  application  of  the  concept  can  assume  a  variety  of  forms  depending 
on  the  type  of  problem  to  be  solved  and  on  the  manner  it  is  modified 
to  account  for  the  constraints.  This  report  is  dealing  with  fixed  and 
constraints  for  some  state  variables  in  general  and  with  the  inventory 
and  production  scheduling  problem  in  particular. 


CHAPTER  2 


THE  METHOD 


1.   A  General  Problem 

An  optimization  problem  in  a  fairly  general  form  can  be  stated  as 
follows:   Determine  0(t)  in  the  interval  tQ  <  t  <  T  SO  as  to  maximize 

*  =  ♦(x(T),  T) 
subject  to  the  constraints 

*  =  *(x(T),  T)  =  0 

3f-  f(x(t),  Jit),  t) 


(1) 

(2) 
(3) 


vith  tQ   and  x(t  )  being  given 
where 


e(t) 


an  mxl  matrix  of  control  variables 


x(t)   = 


f  \M] 


.  *   (t). 


an  nxl  matrix  of  state  variables  which 
results    from  the  choice  of  6_(t)    and  the 


>*P 


an  pxl  matrix  of  terminal  constraint 
functions  each  of  which  is  a  known 
function  of  x(T)  and  T. 


If  an  integral  is  to  be  maximized,  simply  introduce  an  additional 


dx 

~L=   q.(x(t),  6(t),  t)  (U) 

where  q  is  the  integrand  of  the  integral.   Now  x  .,  (T)  is  maximized 

n+1 

with  the  initial  condition  x  ,(t.)  =  0. 
n+l  u 

2.  Concepts 

The  gradient  technique  is  basically  a  method  in  which  the  control 
variable  is  improved  from  a  point  far  away  from  the  optimum  along  the 
gradient  direction.  The  gradient  direction  being  referred  to  is  the 
gradient  that  a  particular  control  variable  has  with  respect  to  the 
given  objective  function. 

Figure  1  is  a  sketch  of  an  optimum  programming  problem.  The  state 
variable  x  must  satisfy  the  functional  relationship 

where  6  is  the  control  variable.  The  functional  relationship  must  be 
satisfied  at  all  points  between  the  two  end  points  t  =  t  and  t  =  T. 

The  problem  of  optimization  basically  arises  when  it  is  necessary 
to  find  the  controls  which,  while  satisfying  the  path,  also  optimizes 
the  objective  function  <Jj.   In  addition  to  optimizing  the  process  and 
satisfying  the  functional  relationship,  it  might  be  required  that  the 
controls  also  meet  certain  additional  final  conditions  on  the  time  variable 
t  or  on  the  final  values  of  the  state  variables. 

To  solve  this  problem  by  the  gradient  technique,  a  certain  sequence 


t  rtQ 


Terminal  Poin't 


l    t=  T 


dx 


X(x,6) 


Pig.l  Symbolic  Sketch  Of  Control  Proble:: 


of  values  is  assumed  for  the  control  variables.  The  gradient  of  the 
objective  function  with  respect  to  the  control  is  determined  at  all  points 
along  the  path  or  trajectory.   According  to  the  gradient  technique,  it 
is  best  to  improve  the  control  variable  along  the  gradient  direction  in 
order  to  reach  the  optimum  as  quickly  as  possible.   The  control  variable 
is  hence  improved  in  the  gradient  direction  by  a  certain  predetermined 
step  size  at  all  points  along  the  path.  At  each  new  point  the  gradient 
is  redetermined  and  a  step  is  again  taken  in  this  new  direction. 
Proceeding  in  this  manner,  an  optimum  is  reached.   After  determining  the 
gradient  direction  and  the1  step  size,  the  new  control  variable  can  be 
computed  from  the  relationship 


e   (t)  =  e  ,.('t)  +  k  II- 

new       old        30  | 


3.   Functional  Equations 

Before  beginning  a  numerical  calculation  based  on  the  gradient 
method,  two  decisions  must  be  made.   First,  a  method  of  calculating  the 
partial  derivatives  making  up  the  gradient  must  be  selected.   Second, 
some  scheme  to  determine  the  step  size  along  the  gradient  must  be 
devised.   The  direction  of  step  along  the  gradient  will  be  determined 
by  whether  a  maximum  or  a  minimum  of  the  objective  function  is  desired. 

Two  main  methods  are  available  to  calculate  partial  derivatives. 
The  simpler  way  is  the  independent  perturbation  of  each  element  of  the 
control  variable.   The  partial  derivatives  can  be  estimated  from  these 
perturbations. 

For  example,  consider  the  process  of  obtaining  the  gradient  with 


respect  to  the  control  0,  the  grid  size  being  A.  The  first  step  is  to 


function  S  .  The  next  step  would  be  to  purturb  9  by  a  small  amount  AS 
and  let  the  new  value  of  the  objective  function  be  S  .  The  gradient  o 
the  objective  function  with  respect  to  the  control  may  now  be  written 


88 

38 


35  I 

vhere  t-         is  the  partial  derivative  of  the  objective  function  with 

36  lt 

respect  to  the  control  variable  at  time  t. 

The  above  procedure  has  the  advantages  of  simplicity  and  of  high 
relative  accuracy  even  in  the  presence  of  nonlinear  effects.  However, 
its  major  disadvantage  is  a  rapid  increase  in  computation  time  as  the 
number  of  time  increments  involved  increases.  For  example,  increasing 
the  number  of  time  increments  from  5  to  10  would  increase  the  computer 
time  required  for  the  calculation  of  the  gradients  by  almost  k   times. 
Thus,  practical  considerations  limit  the  use  of  this  technique  to 
problems  involving  optimization  over  a  relatively  small  number  of 
control  variables  with  a  fairly  small  number  of  grid  points. 

The  second  method  for  calculating  the  partial  derivatives  which 
avoids  the  difficulty  mentioned  above  is  by  setting  up  recurence  relationships 
for  the  calculation  of  the  derivatives. 

Suppose  there  are  n  state  variables  and  one  control  variable. 
Further  assume  that  a  feasible  sequence  of  control  variables  8(t), 
t„  ^  t  <_  T  can  be  obtained  to  start  the  calculations. 


The  performance  equation  as  given  in  Eq.  (3)  may  also  be  written 


x,(t  +  it)  =  x,(t)  +  f.(x.(t),  x„(t),  x  (t),  8(t))  At  (5) 

1  l      1  X     &  n 


where  x.(t)  and  i  =  1,  2,  ...  n  designate  the  value  of  the  i   element 
of  the  state  vector  at  time  t.  Now,  define 


S(x  ,  x  ,  — ,  x  )  =  The  final  value  of  the  objective  function  is 


and  using  the  assumed  control  sequence  9(t), 

Using  Eq.  (5),  then 

S(x, ,  x„,  ...,  x  )  =  S(xn  +  f,  At,  x0  +  fit,  ...,  x  +  f  it).     (6) 
1  c  n      1    X     d         2  nn 

It  is  necessary  to  determine 


as(v  x2,  ...,  xn) 


This  is  the  partial  derivative  of  the  final  value  of  the  objective 
function  with  respect  to  the  control  variable  evaluated  at  time  t.  Using 
Eq.  6,  then 

3S(x,,  x0,  ...,  x  ),     3S(x,  +  f,4,  ...,  x  +  f  A)i 

X       d nj 1    1 n    n  #7\ 

36         |t  38            |t            "I 

How,  by  the  chain  rule 


3(x.  +  f.A) 
1    1 


(6) 


Using  Eq.  (5), 


a(x.  +  f.At)  i     3f. 
l    l  1 


Combining  Eqs.  (8)  and  (9),  gives 


3S 

38 


_   r   db_      i 

.  ~  A  3x.  I   ,   36   . 
t   1=1   l  't+A     't 


(10) 


The  left  hand  side  of  the  above  equation  is  the  desired  partial 

derivative  of  the  objective  function  with  respect  to  the  control  at  a 

3f.  i 


particular  stage  with  a  given  At.   The  partial  derivative 


may  be 


calculated  analytically.   Thus  one  has  a  starting  value  and  a  recurence 
3S 


relationship  for 


However,  before  making  the  calculation,  a 


recurence  relationship  must  also  be  obtained  for  - —     l  =  1,  2, 

3x. 

l    't 

Differentiating  Eq.   6  with  respect  to  x     gives 


3S(xn  ,    >.„,    ...,   x    )    I  9S(x_    +   f.  At,    x„  +  f„At,    ...,    x     +  f  At) 
12                   n            _  1  12  2  n  n 

3XJ  It  3XJ 


(11) 


|    _    y      -      as  | 

it  ~  i=i  3(xi  +  fi  t}   It 


3(x.    +   f.At) 

l  l 

3x, 


(12) 


IE) 


But 


3(x.  +  f .At) 


(13) 


3(x.  +  f.At)  I     3f.  , 

%x   * -  ~    At     i  j  J  (Ik) 


i  =  J  (15) 


Combining  Eq_s.  (11)  through  (15)  results  in 


3S   I     33  ?  3S 


it. 


■  IS-    Tr       "  •     ■     (l6) 

3XJ  [t    3Xj  't+At    i=l  3Xi  't+At  3Xj  It 


The  recurence  relationships  obtained  above  are  essentially  equivalent 
to  those  developed  by  Bryson  and  Kelley  [k ,   6].   However,  the  developments 
due  to  Bryson  and  Kelley  are  considerably  more  complex,  involving 
perturbations  and  adjoint  equations. 

The  derivations  outlined  above  were  based  on  the  situations  where  the 
control  vector  contained  only  one  element.  Problems  with  control  vectors 
containing  several  elements  do  not  cause  any  basic  change  in  the  development. 
However,  there  would  be  individual  recurence  relationships  for  the  partial 
derivative  of  the  objective  function  with  respect  to  each  control  variable. 

Now,  a  relationship  must  be  developed  to  obtain  an  improved  control 


n 


variable  based  upon  the  gradients.  Ao  stated  bcrorc,  the  new  value  of 
the  control  variable  can  be  computed  from  the  relationahip 


i   lt»"8i1>>*1||7  UT1 

new        old  i  't 


Hote  that  e^t)  is  the  value  of  the  ith  element  of  the  control  variable 
at  time  t.   The  scaler,  K,  may  be  thought  of  as  the  step  size  along 
the  gradient.  Thus  the  problem  of  obtaining  a  correction  based  on  the 
gradients  reduces  to  one  of  selecting  a  proper  K.  One  feature  of  K  is 
immediately  obvious:   its  sign  depends  upon  whether  a  maximum  or  a 
minimum  is  desired.  Maximization  problems  require  a  positive  K  and 
minimization  problems  require  a  negative  K. 

A  straight  forward  method  of  obtaining  K  would  be  to  search  over 
all  reasonable  values  and  select  the  one  which  gives  the  maximum 
improvement  in  the  objective  function.   However,  there  is  no  way  to 
specify  the  best  range  of  K  over  which  the  search  should  be  conducted. 
There  is  an  additional  difficulty  of  computation  time.  Since  the 
objective  function  is  defined  at  the  final  time,  evaluation  of  a  trial 
K  requires  a  complete  integration  of  the  performance  equations.  If 
the  process  has  a  large  number  of  state  variables  or  a  large  number  of 
time  increments,  this  integration  would  require  a  great  deal  of  computer 
time. 

An  alternative  to  the  direct  search  for  the  determination  of  K  is 

as  follows.  The  expression  K  rrr~        is  basically  the  difference  between 

1  't 
the  old  and  the  new  values  of  6 .  at  time  t.   In  incremental  form  this 


12 


may  be  written  as 

1 

Suppose  it  is  wished  to  estimate  the  total  change  in  the  objective 

function  due  to  a  series  of  changes  in  8..  One  way  to  obtain  this  estimate 

for  a  process  containing  T  time  increments  would  be  through  the 
approximation 

»♦:!  If-  I   "i(t)  ■  (19) 

t=o  39i  It  x 

Combining  Eqs.  (18)  and  (19)  gives 

-  K  I  flf-  I  f  (20) 


Substituting  the  value  of  K  from  Eq.  (21)  into  Eq.  (IT)  yields 


(21) 


A*  3S 

*  w 

(t)  =  e.   (t)  +  -= ^-L—  .  (22) 

r  old      r      [3S   | 

Jo  i3ei  k 


13 

In  this  equation,  one  must  Delect  a  suitable  value  for  A*|>.   In  other 

voriio ,  it  must  he  decided  how  much  a  change  in  the  objective  function  io 

desired.   If  too  small  a  value  of  A<|>  is  selected,  many  evaluations  of 

the  gradient  will  be  required  to  obtain  the  optimum  while  too  large  a 

value  of  A({>  runs  risk  of  obtaining  no  improvement  at  all.   Notice  that 

only  a  linear  relationship  is  used  in  obtaining  the  gradients.  A  large 

A<(>  can  move  the  new  controls  outside  the  region  where  linearization  is 

valid.  Sometimes  it  may  be  possible  to  obtain  a  scheme  for  adjusting 

A$  during  the  calculations  to  obtain  a  good  balance  between  computation 

time  and  accuracy. 

If  there  are  no  constraints,  one  would  exoect  rr—  to  approach  zero 

l 

as  the  maximum  or  minimum  value  of  S  is  approached.   Since  the  correction 

T   f3s   I  \2 
scheme  requires  division  by  7  \TS~  >  *^e  computation  will  involve 

t=0  L38i  lt  I 

division  by  a  very  small  number  when  the  maximum  or  minimum  is  appraoched. 
The  severity  of  this  difficulty  will  be  discussed  in  later  chapters. 

k.      Functional  Equations  with  Fixed  End  Conditions 

Suppose  that  the  following  condition  must  also  be  satisfied  at  the 
end  point: 

z(xn ,  x„,  x  ,  t)  =  0. 
1       ti       n 

It  is  desired  to  compute  the  influence  of  the  control  variable  6  on  the 
final  value  z . 

The  arguments  used  to  derive  the  relationships  would  be  identical 
to  those  used  in  deriving  the  recurence  relationships.   The  following 
recurence  relationship  for  the  additional  final  condition  can  be  obtained 


14 


°7'  -       V       37'  La  fn->\ 

m  L  "  A  sr  LA  ~o"  L  A  (23) 

't        1=1        1    ' t+A  't 


j    't  j    't+A       i=l  7:  i   't+A      j    't 


with  the  final  condition 


3f.  i 


1?_    -  A5_      |M.  /  M  I  I  -lit  I  fosl 

3xj  It  "  3xj  It  "  ^3t  dt>   It  3xj  It  ' 


Now  let  the  improvement  in  the  control  variable  take  the  form 


i-'t         i  't 


simultaneous  equations: 


^vlfiU^ilif  ItU  ™ 


t=o  ^ae  It  3e  ll 


*■  •  *^"X  f*  I*  *  I;  ]  *  *»  J-  te  L 1*  (23) 


Due  to  the  difficulties  in  finding  the  initial  feasible  trajectory, 
the  term  Az,  which  represents  the  deviation  from  the  desired  final 
condition,  can  be  thought  of  as  a  correction  in  the  final  value  of  the 


15 


subsidiary  condition. 

It  is  to  be  noted  that  the  expression 


2  36.  L 
l  't 


in  Eq.  (26)  is  a  penalty  function  imposed  on  the  improvement  of  the 
control  variable  due  to  the  auxiliary  final  condition  which  must  be 
satisfied  by  the  optimal  sequence  of  the  control  vector.  The  penalty 
function,  in  general,  reduces  the  rate  of  approach  to  the  optimum. 
Summarizing  the  above  discussion,  it  is  seen  that  if  some  end 
conditions  at  the  final  time  have  to  be  satisfied,  Eq.  (17)  must  be 
modified  to 


6i        <t>-6  <t)+     Klff-         +K2ff-  (29) 

new  old  l  l    't  l    't  ' 

-    I 

where  the  constants  Kn  and  K„  and  the  expression  for  -rr-—         are  computed 
1      2  36.  |t 

from  the  relationships  specified  in  Eqs.  (23)  through  (28). 


5.   Computational  Procedure 

Knowing  the  recurence  relationships,  the  step  to  step  computational 
procedure  for  a  problem  without  given  end  conditions  are: 

1.  Assume  a  sequence  of  control  variables. 

2.  Using  the  performance  Eq.  (5),  determine  the  sequence  of  state 

variables. 

3f.  I       3f.  | 

3.  Evaluate  t——    and  - —    for  every  t  using  the  numerical 

36   |t      3Xj  |t 

values  of  the  state  vectors  and  controls. 
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k.      From  the  final  values  of  the  state  vector  and  the  controls, 

evaluate  —' —  at  the  end  of  the  process. 

l 

5.  Calculate  rr—    by  backward  recursion  of  Eq.  (16). 

86j  it 

6.  Calculate  —    from  Eq.  (10). 

36  |t 

T.   Calculate  the  new  control  from  Eq.  (22). 

8.  Repeat  Steps  2  through  7  until  the  gradient  is  so  small  that 
further  improvement  is  not  significant. 


In  case  the  problem  involves  satisfying  some  end  conditions  on  the 
state  variables,  the  only  difference  in  the  computational  procedure 
would  be  in  Step  7  where  the  improved  control  would  now  be  computed  from 
Eq.  (26)  instead  of  Eq.  (22). 


17 


CHAPTER  3 

APPLICATIONS 

Three  problems  relating  to  the  optimization  of  management  systems 
have  been  solved  by  the  gradient  technique  in  this  report.  The  first  is 
a  simple  problem  in  the  field  of  production  control.  Here  both  the 
initial  and  the  final  inventories  are  given.  Thus  this  is  an  optimization 
problem  with  one  fixed  end  condition.  The  sales  in  this  case  are  assumed 
to  be  known.  The  second  problem  considers  the  diffusion  model  of 
advertising.  The  problem  is  to  find  the  optimum  advertising  for  the 
maximum  profit  with  a  given  production  rate.  The  model  also  controls 
the  inventory  level.  The  final  problem  considers  both  production  and 
sales  as  variables  with  the  operating  temperature  controlling  the 
production  rate  and  the  diffusion  model  being  used  for  the  advertising. 

3.1  A  Production  and  Inventory  Control  Problem 

3.1.1  The  Model  "  ■  .  • 

Consider  the  solution  of  a  simple  problem  in  the  field  of  production 
and  inventory  control  using  the  gradient  technique.  Further  consider  the 
sales  rate  Q(t)  to  be  known  with  certainty  and  that  the  rate  of  change 
of  the  inventory  level  I(t)  is  given  by 


^-=p(t)-Q(t)  (30) 


where  p(t)  is  the  production  rate  at  time  t.  The  problem  is  to  minimize 
the  cost  function 


IB 


CT  =  ^   k(Kt))2  +  C  exp  (p  -  p(t))2]dt  (31) 

0   It  r        m  j 

where  CT  is  the  total  cost  of  inventory  and  production.  C  is  the 
inventory  carrying  cost.   C   is  the  minimum  production  cost  which  occurs 
when  the  production  rate  equals  p  which  can  be  considered  as  the 
capacity  of  the  manufacturing  plant.  Since  the  plant  is  designed  for 
capacity  p^,  an  increase  in  capacity  may  require  additional  equipment 
and  manpower  and  thus  can  be  very  expensive.  On  the  other  hand,  a 
decrease  in  capacity  can  be  equally  expensive  because  of  maintenance 
of  unused  equipment  and  because  of  manpower  which,  due  to  contract 
agreements,  cannot  be  easily  reduced. 

Sales,  a  known  quantity,  is  given  by  the  linear  relation 

Q(t)  =  a  +  b(t)  .  (32) 

The  initial  inventory  is 

1(0)  =  C  (33) 

It  is  desired  that  the  inventory  level  at  final  time  t  =  1  be  I(T)  =  D. 

Reformulating  this  problem  according  to  the  notations  used  previously 
results  in  x^t)  =  l(t)  and  6(t)  =  p(t).  Equation  (|6)  now  becomes 

dxx(t) 

— jr— «  6(t)  -  a  -  b®  (3U) 


xl(0)  =  C  (35) 


xx(T)  =  D.  (36) 
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x2(t)  =  /*  [CI(x1(t))2  +  C  exp  (pm   -  p(t))2]  4t.  (37) 


x2(T)  =  CT 


-jj|-  CI(x1(t))2  +  C  exp  (pm  -  p(t))2  (38) 


x2(0)  =  0.  (39) 

Equations  (3^)  and  (38)  are  the  desired  differential  equations 
corresponding  to  Eq.  (3)  with  n  =  2.   The  initial  conditions  are  given 
by  Eqs.  (35)  and  (39).   The  function  to  be  minimized  is 

*  =  x2  .  (U0) 

The  terminal  conditions  are 

<|i  =  t  -  T  =  0  (1*1) 

and 

z  =  x  (T)  -  D  =  0  .  .  (1(2) 

3.1.2  The  Recurence  Equations 

Considering  t  as  the  third  state  variable,  the  variational  equations 
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can  bo  obtained  easily.   From  Eq.  (10), 


f§    =  ~      A  -  2  |f-      [C  (p  -  6}  exp  (p  -  6)2]  (U3) 

3*  U   axi  l-t+a    3x2  lt+4  p  ffi        m      U 


Using  Eq.  (l6)  yields 


is  I       is    :       ^      is    I       „  „, , 

3X1    >t        3X1    U+fi  3X2    U+A     J      X    U 


(U5) 


The  terminal  conditions  are  obtained  from  Eq.  (9)  give 


||-   -  0  (1,6) 

3X1  'T 


ff-   -  1  -  (1*7) 

2  't 

An  equation  for  rr  can  be  obtained  from  Eq.  (8).   However,  for  the 

present  problem  this  equation  is  not  needed. 

The  variational  equations  for  the  second  constraint 

z  =  x  (T)  -  D  =  0 

can  be  obtained  in  a  similar  manner.   From  Eq.  (23), 
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From  Eq.    {2k) 


3xi  It      3xi  U+a         3x2  U+a    *    x|t 


lr     -fx~       •  (50) 

3X2    It        3X2    U+A 


The  terminal  conditions  can  be  obtained  from  Eq.  (25] 


3z    _  .  ,  , 

3x"~    =  2  (5l' 

Xl  <T 


3xl  It 

From  Eqs.  (Uo)  and  (1*2),  then 


2  U   3X2  It+A 


(53) 


Equations  (39).  (kl)   and  (1*3)  give 


3z     _  3z_  I  ,   _  .  ,,,  , 

3X1  U   3X1  't+A 


Substituting  the  values  obtained  in  Eqs.  C*3)  through  (51)  into  Eqs.  (27) 
and  (25),  the  values  of  the  constants  K  and  K  can  be  calculated.  By- 
substituting  these  values  of  K  and  K  into  Eq.  (29)  the  improvement  in 
Az  is  obtained.  The  value  of  Az  is  (x 
improvement  in  the  objective  function. 
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3-1.3  Numerical  Results 

The  numerical  values  used  are 

a  =  2  D  =  9.25 

*>■*  1  C  =  0.001 

p 
C  =  5  pm  =  5 

Cj  =  0.1  T  =  1 

A  =  0.01 

This  problem  was  solved  on  an  IBM  360-50  computer.  The  convergence 
rate  of  the  control  variable,  the  production  rate,  is  shown  in  Fig.  (2). 
The  convergence  rates  of  the  inventory  level  and  the  cost  function  are 
shown  in  Figs.  (3)  and  (It),  respectively.  The  Runge-Kutta  integration 
formula  was  used  to  integrate  Eqs.  (39)  and  (38).   The  step  size  used 
was  0.01,  which  is  the  same  as  the  A  value  used.  A  value  of  A(J>  equal 
to  -  0.1  was  used  for  the  first  25  iterations  and  values  of  A<*>  =  -  0.01 
and  A<)>  =  -  0.001  were  used  for  26  to  72  and  73  to  126,  respectively. 

The  last  part  of  the  production  rate  curve  is  very  insensitive 
to  the  cost  function  x2-  Only  five  iterations  were  required  to  get  a 
cost  very  near  the  optimum.  However,  the  curve  for  the  production  rate 
at  the  fifth  iteration  is  far  from  the  optimum  one.   (See  Fig.  2).  The 
cost  C^  at  the  fifth  iteration  was  5.25;  the  minimum  cost  was  5.17,  a 
decrease  of  only  l.lltjS. 

The  convergence  rate  from  the  fifth  iteration  to  the  optimum  is 
very  slow.  Approximately  90  iterations  were  required  to  improve  the  cost 
from  5.25  to  5.17.  This  difficulty  cones  from  the  fact  that  the  gradient 
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• 
is  very  small  near  the  optimum.  The  end  condition  required  of  the  inventory 
vas  satisfied  at  the  second  iteration,  and  the  inventory  converged  to  the 
required  final  value. 

3.2  An  inventory  Control  Problem  with  Advertising 

3.2.1  The  Model 

Consider  an  inventory  control  model  as  shown  in  Fig.  5-  With  the 
production  rate  given,  the  problem  is  to  balance  the  cost  of  advertising 
and  the  inventory  level  for  maximum  profit.  The  variables  involved  are 
production,  inventory,  sales  and  advertising.  Production  can  be  considered 
a  function  of  time.  Let  production  at  time  t  be 

■  P(t)  =  A  +  bt  (1*6) 

where  A  and  b  are  constants. 

The  rate  of  change  of  inventory  I(t)  is  the  difference  between  the 
production  and  sales  at  time  t.  If  Q,(t)  represents  the  number  of 
customers  at  time  t  and  C  represents  the  number  of  times  bought  by 
each  customer,  the  rate  of  change  of  inventory  level  may  be  represented 
by 

^-=p(t)  -CqQ(t)  (VT) 

To  determine  the  rate  of  change  of  customers  (informed  persons) 
at  time  t,  a  diffusion  model  incorporating  advertising  will  be  used. 
Consider  a  group  of  people  in  which  only  a  certain  number  possess  a 
particular  piece  of  information.   Suppose  that  the  total  number  of  persons 
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in  the  group  under  consideration  remains  constant  and  that  diffusion  of 
information  only  occurs  through  personal  contact.  The  number  of  contacts 
made  by  an  average  person  in  an  arbitrary  unit  of  time  is  given  by  a 
contact  coefficient.  In  a  contact,  the  contactee  receives  information 
if  he  does  not  have  it;  if  he  already  has  it;  the  contact  is  wasted  so 
far  as  increasing  the  number  of  people  who  have  the  information  is 
concerned. 
Let 

Q(0)  =  the  number  of  informed  participants  in  the  group  at 
time  0. 
N  =  the  total  number  of  participants  in  the  group. 
Q(t)  =  the  number  of  informed  participants  at  time  t. 

am 


1  -  "    ■  proportion  of  uninformed  persons  in  the  group  at  time  t. 

C  Q(t)dt  =  contacts  made  during  a  time  interval  dt. 
1 

The  increase  in  the  total  number  of  informed  persons  during  a  short 
interval  At  is  obtained  by  multiplying  the  number  of  contacts  by  the 
proportion  of  persons  who  do  not  possess  the  information,  since  only 
contacts  with  uninformed  persons  leads  to  an  increase  in  the  informed 
members . 

dQ(t)  =  CQ(t)[l  -  ^  ]At  (W) 

The  differential  equations  is 


2J 


«t4l.fl,(t)U-^l]  (u9) 


Suppose  that  the  information  in  the  model  given  above  is  about  the 
product  of  a  firm  and  assume  that  the  firm  can  influence  the  number  of 
contacts  by  spending  money  for  advertising.  In  particular,  it  can  in- 
crease the  number  of  contacts  made  by  the  informed  persons  (above  the 
ones  included  in  C)  by  an  additional  number  per  unit  of  time.  Then 

aQ(t)dt  =  number  of  additional  contacts  made  in  the  time  interval  dt . 

Hence  the  differential  equation  becomes 

«J*1.  [C  +  a]Q(t)  [l-^1]    .  (50) 

The.net  profit  can  be  obtained  from  the  equation: 

Net  Profit  =  Revenue  -  Cost  of  holding  inventory  -  Cost  of 
advertisement. 


Since  C  Q(t)  units  are  sold  at  time  t,  the  revenue  is  FC  Q(t).   If  I 

q  q  ir 


p 

cost,  the  cost  due  to  inventory  is  CT[I  -  I(t)]  .  The  cost  of  advertising 

1  m 
p 
is  C  Q(t)  A  (t)  and  the  total  net  profit  over  the  duration  of  the  process 

is 


J  =  /   [FC  Q(t)  -  C  [I   -  I(t)]2  -  C  Q(t)a2(t)]dt.  (51) 

q     1  1   ">  A 

3.2.2  The  Recurence  Equations  ' 

The  state  variables  are  l(t)  and  Q(t)  and  the  control  variable  is 
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the  additional  contact  coefficient  a(t).  The  objective  is  to  maximize 
the  total  net  profit  J.  The  differential  equations  representing  the 
process  are  Eqs.  (1*6),  (I47)  and  (kg)   with  the  profit  represented  by  J. 

To  reformulate  the  problem  in  terms  of  the  derivations,  let  l(t)=x  (t), 
Q(t)  ■  x2(t)  and  a(t)  =  6(t).  In  addition  the  following  variable  can  be 
introduced: 

,"*  ?  ? 

x,(t)  =  /.  [FC  Q(t)  -  CT[I  -  I(t)r  -  CA  Q(t)a*(t)]dt  (52) 

i  1    q       1  m  A 


-rr1*  [FC  Q(t)  -  CUl  -  I(t)]2  -  C  Q(t)a2(t)]  (53) 

ut      q       1  n  A 


with 


and 


x3(0)  =  0 


x3(tf)  =  J  . 


The  objective  now  became  the  maximization  of  <{>  =  x  (t_)  with  the  terminal 
condition 

*  =  t  -  T  =  0  . 
Using  Eq.  (12),  these  recursive  relationships  are  obtained: 


3S_  I   _  3S_  I    +  fas_  I    !£l|+3S_|    af2  I     3S   I     3f3  I  ) 

3xl  It  "  3X1  If*   K  lt+A  3xl  It   3x2  If  4  3X1  It  +  3X3  lt+A  3X1  I J  ^ 


(5M 
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^  t"^LMid  <2CI(I»-X1»I 

1    't  1    't+A        "■     3     t+A 


3S      I      _   3S_  I  +    I" 3S_  I  ^1.1      +3S_!  !I2|+1S_|  !!l   I      |    A 

8X2    It  "   3X2    U+A        K    I t+A   5X2    It        3X2    U+A   3X2    U        3X3    't+A   3X2    U  > 


(56) 


I- 1     ♦lr|     (c+eja-^  +  ff-l    im/)1. 

3X1    I  t+A        3X2    U+A  N  3X3    't+A  A       > 


(57) 


and 


3S       3S       ^   3S 


3x3  It   3x3  I t+A    l8xl  I t+A  3x3  It   3X2  I   A  3x3  It 


3S   I     3f3 


3x3  't+A  3x3  U 


3S   [      3S_  | 
3x3  U+A  =  3X3  U  ' 


(59) 


Terminal  conditions  obtained  from  Eq.  (16)  are 


1—    -  0  (60) 

3xi  It 


ft-."0  (6l) 

2  >T 
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3x3  |T  "  X>  (62) 


To  find  the  recursive  relationship  for  the  gradient  of  the  objective 
function  with  respect  to  the  control  variable,  there  is  the  general 
relationship 

as  I   K1   as  I   3fi  I 

39  L   A  3x.  L.  ~Te   L  (63) 

't   1=1   i  't+A     't 

2 

(3S   I  x?    3°  1 

Equations  (ll)  through  (21)  are  the  required  recursive  relationships 
for  the  solution  of  the  above  problem  by  the  gradient  technique. 

3.2.3  Computational  Procedure 

The  numerical  values  used  were 

CJ   =  0.15  A  =  70 

CA  =1.5  b  =  100 

F  =  10.0  C  =  2.0 

B  =  150  C  =  1.0 

q 

Im  =  50.0  At  =  .01,  T  =  1 

The  values  of  A*  used  were:   A$  =  1+0  for  first  17  iterations,  and  A*  =  .05 
for  the  remaining  iterations.  With  the  initial  conditions  as 


X  (0)  =  20.0,   x  (0)  =  20.0   and   x  (0)  =  0, 

the  step-by-step  procedure  followed  for  obtaining  the  solution  was: 

1.  Integrate  Eqs.  (1)7),  (50)  and  (51). 

2.  Obtain  the  end  conditions  for  —  ,  gj-  ,  J^~  ,     from  Eqs.  (60) 

through  (62). 

3.  Calculate  the  values  of  ||—  ,  ~—  and.  |^—  at  the  101  grid 

points  by  means  of  backward  recursion  of  Eqs.  (55)  through  (59:. 

9S 
It.  Calculate  the  gradient  of  the  control  variable  j^  at  the  101 

grid  points  by  means  of  backward  recursion  of  Eq.  (6U). 

5.  Calculate  the  improvement  in  the  control  variable  8  by  the 
relationship 

3S/36 
new    old      T  ,~  2 

y   (— ) 

6.  Repeat  Steps  (l)  through  (5)  till  no  further  improvement  could 
be  made. 

3. 2. It  Discussion  of  Results 

Using  an  initial  guess  of  2.5  for  the  advertising,  the  control  variable 
converged  to  the  optimal  in  80  iterations.  The  convergence  rate  of 
advertising  is  shown  in  Fig..  6.  The  profit  function  had  a  value  of  25-0 
with  the  assumed  controls  and  converged  to  the  optimal  value  of  530.0  in 
80  iterations.   It  is  seen  that  the  rate  of  convergence  during  the 
initial  stages  of  iteration  was  much  faster  than  during  the  final  stages. 
The  reason  for  this  is  that  the  gradient  is  small  and  hence  improvement 
is  also  small  at  the  final  stages  of  iteration.  The  trajectories  of  the 
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state  and  control  variables  are  given  in  Figs.  6  through  9. 

3. 3. A  Production,  Inventory  Control  and  Advertising  Problem 

3.3.1  The  Model 

Consider  the  manufacturing  process  shown  in  Fig.  10.  The  raw 
material  is  fed  into  two  reactors  in  series  and  the  Arrhenius  reaction 
rate  expression 

k  =  G  exp  (-  g|) 

will  be  used  for  the  rate  constants.  The  reactions  are  first  order  and 
can  be  expressed  as 

A  $  B  &  C 
where  k  is  the  reaction  rate  constant,  G  is  the  frequency  factor  constant, 
E  is  the  activation  energy  of  the  reaction,  R  is  the  gas  constant  and  T 
is  the  temperature.  Material  B  is  the  desired  product  for  which  inventory 
and  advertising  are  assumed. 

The  transformation  equations  for  the  two  reactions  are 


i  dt  =  q(xo  "  Xl'  "  Vl  V  Xl  (65) 


vi  ir  ■  *(y0  -  yi>  -  vi  \  yi +  vi  \  vi  (66> 


V2  dT  ■  *(X1  -  X2>  "  V2  \   X2  (6T) 


V2  ^  =  q(y,  -  y2)  -  V2  K     y2  ♦  ^   k   x,,  (68) 
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respectively,  q  is  the  flow  rate,  k   and  k   represent  the  reaction  rate 

i      i 

for  the  first  and  second  reactions,  respectively,  and  x  and  y  are  the 
concentrations  of  A  and  B,  respectively. 

If  C  represents  the  number  of  items  bought  by  each  informed  person, 
the  change  in  inventory  l(t)  can  be  represented  by  the  differential 
equation 

«g*l  =  qy2  -  Cq  K(t)  (69) 

where  K(t)  is  the  number  of  informed  persons  at  time  t. 

To  determine  the  sales,  the  diffusion  model  of  advertising  discussed 
in  Section  3-2  will  be  used  and  the  differential  equation  is 

SfJ^-  [C  ♦  a(t)]  K(t)  [1  -  *]&•]  (70) 

The  profit  equation  can  be  written  as 

Net  Profit  =  Income  from  sales  of  A  +  Income  from  sales  of  B 
+  Income  from  sales  of  C  -  Cost  of  Inventory 
-  Cost  of  Advertising  -  Cost  of  Production 


the  profit  equation  can  be  written  as 


J  =  /   [C^KU)  +  C2qx2  +  C3q(l-x2-y2)  -  Cjd^-Kt )  f  -   CA[a(t)K(t)! 
0 

"  CTtTlm"Tl)2  +  fri-*2>2n*t- 
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where  C  ,  C  ,  C  ,  C  ,  C.  and  C  are  the  per  unit  costs  of  B,  A,  C, 

inventory,  advertising,  and  production,  respectively. 

Introducing  an  additional  state  variable  x  now  gives 

,*  ?  ? 

Xj(t)   =   /    [C1C<iK(t)   +  C2qx£(t)    +  C3q(l-x2-y2)    -  C^I^T  -  Cje^r 

-  cT[(Tlm-T1)2  +    (T1-T2)2]]dt  (71) 


with  x  (0)  =  0  and  x  (t  )  =  J.  Representing  K(t)  by  x, (t)  and  I(t)  by 

x,(t),  a  problem  involving  six  state  variables,  x  ,  y  ,  x  ,  y  ,  x  and  x, 

and  three  control  variables,  T  ,  T  and  a(t),  is  involved.  Let  V  =  V  =  V; 

now  the  differential  equations  representing  the  process  may  be  summarized 
as: 

dx^   q(xQ  -  xx) 

dt  "     V      "  Ka1  Xl 

dyl   q(y0  -  yl> 
dt       V 

dX2   q'Xl  ~  ~A2I 


dy2   q(yx  -  y2) 


dx3 

tt~  =  ^2  "  Cq  Xk 
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.  2 

dx,  x. 


dx 


V^/  +  (Tl-T2)2] 


vith  the  given  initial  conditions  x  (0)  =  x  ,  y.(0)  =  y  ,  x  (0)  =  x  , 

y2(0)  =  y°,  x3(0)  =  x°,  x^to)  =  x°  and  x5(0)  =  0. 

The  problem  is  to  maximize  $  =  x  (T)  subject  to  the  end  condition 
constraint  iji  ■  t  -  T  =  0. 

3.3.2  The  Recurence  Equations 


For  the  end  condition  of  the  slope  of  the  objective  function  with 
respect  to  the  state  variables,  this  equation  exists: 


!§_  I   -  lt_  I   _  fit  /  Jit)  I   14-  |      1  =  i2       7 

3xj  If      3xj  It      l"'  atJ  It  3xj   It         J      1,2>  •••>T- 


For  the  present  case,  then 


3S 
3x~~        =  °  (T2) 

1    'T 


w:lm0  (73) 


2    'T 


"  °  ilk) 


§-    =°  (75) 


35  '  =0  (76) 


3x3  't 


3xit  't 


||-   -  1.  (76) 

3x5  't 


The  recursive  relationship  as  derived  in  Chapter  II  are: 


dS       _  dS   [      _   r   dS         1      . 

3x,    ~  3x,  I  •,    .*,  3x.   ...  3x,   , 

J  't     J  't+A   1=1   l  't+A   j  't 


Using  this  equation  for  the  state  variables,  then 


K  l«M    Vl    "l   U 


3yl    I  t+A        al    't        3X2    U+A  2 


1         \      U 


|f-  (q/Vj|    A 

'2    't+A 


3S  3S  [3S     I  ,  .  I 

sr     =  ir  \     +  tar       (_  q/v?  -  k«  n 

dX2    U        3X2    U+A        ^3X2    U+A  2  a2    I 


3S  3S      I  3S      I  ,         ,  ,| 

*y2    U        *y2    U+A  3x5    U+A  2         ^2    U 


as    |    _  3s_ 


3    U+A      <*        3xl,    U+A  3 


The  recurence  relation  for  the  control  variables  are: 


.„   I     n+1 
3S      _   y   c 

as      L 


it. 


36,  !     ,',  3x.  I   .  36,  , 
J  't   1=1   l  't+A   j  't 


Applying  this  equation  to  Eqs.    (l)   through   (7)   gives 
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"  Ua    »        +   37"  ^-0  4)1    A 

2    U+A        a2    U        dx5    U+A        d  6 


(81) 


3x 


1  +  i7"  «-C31»        A 

3   't+A  dx5   U+A  J 


(82) 


3S 


„   3S_   J  +    (3S_   | 


(2CT(I      -  xj)  A 


t_3x3U+A  Tl3x5lt+i-I-»-X3"|t 


(83) 


(1"^)+^L+A(c^-2Ca6Mi 


(81.) 


3S  _   3S_ 

3xr    I      ™   3x^ 

5    't  5    't+A 


(85) 
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k 
3S_  I  fas      I  ,  _^l.|  3S_    I 

3T1    It        K    If*  1T1      It        8yl    If  A 


3k  3t_ 

(    _!i  _A5 1   +  ii_  I 

lXl     JT.  "  yl     3T/    .        3xB    L. 

1  1   't           5     t+A 


(-  CTC-  2   (Tlm  -  Tx)       +  2(TX  -  T2)      ))|    A  (86) 


3k  3k 

a2, 1    3S   I     f   _! 
3T2  U"  i3x2  If 4  %"  "23T2   U    *2  If  A   2   3T2 


•3S        3S   I     ,        2  J     3S   I     ,      2 


^2,  I  .  as 


and 


-y2lf)    +ff-      (-  C  (-  2(T  -T  )   ))   A        (87) 
2  3T2   U   3X5  't+A     T      X  2  't   ' 


2 


The  improvement  in  the  control  is  given  by 


a*   3S 

e    (t)  _  8    (t)  ♦  — J— L-         J  =  1,  2,  3  (89) 

Jnew       Jold      c  3S_  I 

t£o  38j  it 


where  A<(>   is  the  desired  improvement  in  the  objective  function  due  to  the 
jth  control. 
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Equations  (65)  through  (89)  are  the  desired  equations  for  the  solution 
of  the  problem  by  the  gradient  technique.  ."It  should  be  noted  that  since 
the  decision  vector  is  multidimenstional,  individual  improvements  A$  are 
suggested  for  each  control. 

3.3.   Numerical  Results 

Based  on  the  model  stated  above,  the  problem  was  solved  using  different 
parameters  and  different  starting  values.   The  parameters  and  starting 
values  used  are  summarized  in  Table  1.  The  initial  conditions  for  the 
seven  state  variables  used  are  shown  in  Table  2.  The  results  are  discussed 
in  the  following  sections  for  each  of  the  problems  shown  in  Table  1. 

Problem  la:  It  can  be  seen  from  Fig.  11  that  the  optimal  temperature 
profile  for  T  had  a  value  of  about  362°K  at  time  t  and  339°K  at  tf . 
The  concentration  of  A  in  Fig.  12  fell  to  a  value  of  .1*58  at  time  O.65 
and  raised  to  0 . U8U  at  t  .  The  concentration  of  B  in  Fig.  13  arises  to 
0.1*65  at  t=0.65  and  falls  to  0.1*5l*  at  t  .  The  profit  as  can  be  seen 
from  Fig.  20,  with  the  initial  controls  was  $T2.0l»  and  at  the  20th  iteration 
it  was  to  $97.05,  an  increase  of  3l*.5#.  However,  after  the  20th  iteration 
it  took  another  200  iterations  for  the  profit  to  reach  its  optimal  value 
of  $107.01*.  This  slow  increase  of  only  l.OJf  in  200  iterations  was  due 

to  the  slow  convergence  rate  near  the  optimal.  At  the  20th  iteration 

3S  —  3  — S 

the  sum  of  (-r=-)  was  0.1*  x  10  .  This  sum  was  0.2  x  10   for  the  optimal 
3T1 

profile.   Since  the  value  of  the  gradient  is  very  small  at  this  point, 


3S  2  —1* 

sum  of  (r^— )   at  the  initial  iteration  was  0.3  x  10   .   At  the  20th  it 
3T 

e       T  -3  -1* 

was  0.6  x  10     and  at  the  optimal  it  was  0.3  x  10   .   It  can  be 

seen  in  this  case  that  the  gradient  at  the 


Table  1.   Parameters  and  Initial  Approximations 
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INITIAL 
GUESS 


PARAMETERS 


R 
1 

Vl 

Ga 

E 
a 

V2 

\ 

C 

q 
c 

N 
Cn 


Prob.   la 

Prob.   lb 

Prob.   2 

2.0 

2.0 

2.0 

60.0 

60.0 

60.0 

12.0 

12.0 

12.0 

•535X1011 

.535xl0U 

.535xl0n 

18000.0 

18000.0 

18000.0 

12.0 

12.0 

12.0 

,lt6lxl018 

.Wixio18 

.l*6lxl018 

30000.0 

30000.0 

30000.0 

1.0 

1.0 

1.0 

,     1.0 

1.0 

1.0 

100.0 

100.0 

100.0 

5.0 

5.0 

5.0 

0.0 

0.0 

0.0 

0.0 

0.0 

•  0.0 

1.0 

1.0 

1.0 

0.0002 

0.0002 

0.01 

0.005 

0.005 

0.005 

0.53 

0.53 

0.53 

0.1*3 

0.1*3 

0.1(3 

3l»0.0°K 

31*0. 0°K 

31*0. 0°K 

330. 0°K 

3U5. 0°K 

31*5. 0°K 

330. 0°K 

31*5. 0°K 

31*5. 0°K 

8.0 

10.0 

3.0 

Table  2.  Initial  Conditions  for  the  State  Variables 
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Variable 

Prob.  la 

Prob.  lb 

Prob.  2 

x1(0) 

0.53 

0.53 

0.53 

yx(o) 

0.1*3 

0.1(3 

0.U3 

x2(0) 

0.53 

0.53 

0.53 

y2(o) 

0.1*3 

0.1(3 

0.1(3 

x3(o) 

8.00 

8.00 

8.00 

Xfc(O) 

0.10 

0.10 

1.00 

x  (0) 

0.00 

0.00 

0.00 
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20th  iteration  was  steeper  than  the  initial  gradient  and  the  gradient 
at  the  optimal  vas  not  very  different  from  the  gradient  at  the  initial 
guessed  control.  Since  the  problem  vas  to  optimize  the  process  with 
respect  to  the  three  controls,  it  was  not  necessary  to  have  the  same 
gradient  for  all  the  controls.  The  problem  was  to  adjust  the  step  size 
for  each  control  so  that  the  optimum  could  be  obtained  with  the  minimum 
amount  of  computation. 

In  the  case  of  the  advertising,  (please  refer  to  Fig.  IT)  the 
gradient  at  the  initial  iteration  was  3.0.   This  indicates  that  the 
initial  guess  was  fairly  far  from  the  optimum  and  is  also  evident  from 
the  plot  in  Fig.  17.  At  the  20th  iteration  the  gradient  was  0.2;  at  the 
optimal  it  vas  0.3  x  10  .  From  Fig.  18  it  is  seen  that  the  major 
change  in  the  number  of  informed  persons  occurred  during  the  initial 
stages  of  the  process.  This  fact  is  also  evident  from  Fig.  17  where 
the  additional  contact  coefficient  rises  to  16.95  at  t  =  0.U  and  then 
falls  to  a  value  of  0  at  the  final  time.  It  should  be  noted  that  a  non- 
negativity  constraint  was  used  for  the  additional  contact  coefficient. 
The  plots  of  temperature  and  concentrations  of  products  A  and  B  versus 
time  in  the  second  reactor  can  be  seen  from  Fig.  ih   through  16. 

Problem  lb:   In  this  problem,  the  same  parameters  were  used  as  in  problem 
la  except  that  different  starting  values  were  used  for  the  controls. 

From  Fig.  30  it  is  seen  that  the  profit  at  the  initial  iteration 
vas  $55-90  and  at  the  20th  iteration  it  increased  to  $80.50,  an  increase 
of  U5i«.   However,  after  the  20th  iteration  it  required  180  further 
iterations  to  improve  the  profit  to  $107. 16,  which  is  the  same  result 
as  in  problem  la.  The  percentage  improvement  vas  only  2©5  in  180 
iterations. 
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The  concentration  of  B  in  the  second  reactor  as  can  he  seen  from 
Fig.  26  is  seen  to  converge  quite  rapidly  to  a  point  dorse  to  the  optimal 
However,  the  rate  of  convergence  from  this  point  to  the  optimal  is  very 
slow. 

For  the  temperature  in  the  first  reactor,  the  value  of  the  sum  of 

V3T~ 


? 
3S  — "^  -5 

(t=— )   at  the  initial  iteration  was  0.9  x  10   ;  this  value  was  0.2  x  10 


1 
at  the  optimum.  A  value  of  A<t>  equal  to  0.1  was  used  for  the  first  fifty 

iterations  and  values  of  A<t>  =  0.01  and  0$  «  0.001  were  used  for  51  to 

99  and  100  to  150,  respectively.  Thereafter  the  value  of  A$  was 

successively  reduced  and  the  optimal  was  reached  with  a  A<J  value  of 

0.001  for  T  .  As  has  teen  stated  previously,  the  initial  guess  for 

advertising  was  far  from  the  optimal.  This  necessitated  the  use  of  A<i>  = 

1.0  for  the  first  70  iterations;  thereafter  A<J  =  0.5  was  used  for 

iterations  71  -  100  and  Aiji  =  0.01  was  used  for  iterations  101  through 

150.  The  value  of  A*  for  the  temperature  in  the  second  reactor  was 

0.1  for  the  first  fifty  iterations  and  an  optimal  was  reached  with 

A*  =  0.001. 

The  factor  that  governs  the  choice  of  a  proper  value  for  A<J  is  the 

relative  position  of  the  current  control  with  respect  to  the  optimal 

control.  In  nearly  all  practical  situations,  the  proper  choice  of  A$ 

must  be  obtained  by  a  trial  and  error  procedure.  However,  if  a  certain 

extimate  for  the  location  of  the  optimal  control  exists,  the  required 

computation  to  obtain  the  optimal  could  be  greatly  reduced.   The 

temperature  and  concentration  profiles  in  the  two  reactors  are  shown  in 

Figures  21  through  29. 
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Problem  2 

The  model  used  for  this  problem  was  the  same  as  that  used  for  problems 
la  and  lb  except  for  different  values  of  the  parameters  and  initial  guesses 
for  controls. 

The  profit  with  the  initially  guessed  control  was  $U.10  and  the 
optimal  profit  was  $66.05.  Again  it  can  be  seen  that  the  convergence 
rate  is  very  slow  when  the  current  result  is  near  the  optimum.  The  profit 
at  the  60th  iteration  was  $27.50.  Another  200  iterations  were  required 
to  reach  the  optimal.  The  temperature,  concentration,  inventory,  sales 
and  profit  profiles  are  plotted  in  Figures  31  through  1*0. 
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CHAPTER  U 
CONCLUSION  AND  DISCUSSION 

Although  the  three  problems  solved  in  Section  3.3  have  different 
parameters  and  different  starting  values,  they  have  certain  characteristics 
in  common.  To  study  the  common  characteristics,  consider  the  detailed 
behavior  of  the  system  at  any  particular  iteration. 

At  time  t  ,  the  number  of  informed  persons  in  the  group  is  lov  so 
the  number  of  items  sold  is  low.  There  is  a  large  number  of  uninformed 
persons  and  hence  advertising  is  high.  Since  the  system  is  already 
in  production,  there  is  a  tendency  for  stock  to  go  into  inventory.  As 
time  goes  on,  the  sales  or  the  number  of  informed  persons  increases  and 
the  additional  contact  coefficient  also  increases  because  it  is  profitable 
to  do  so.  The  production  rate  of  B  rises  because  it  has  to  meet  the  sales 
requirement  as  well  as  to  maintain  the  required  inventory  level. 

There  comes  a  stage  when  the  sales  is  much  more  profitable  than  slight 
variations  in  the  inventory  level  and  the  inventory  begins  to  fall.   In 
this  case  this  comes  at  about  t  =  0.55.  However,  when  the  inventory  has 
become  sufficiently  low,  the  inventory  cost  in  the  total  profit  equation 
becomes  important.  It  should  be  noted  that  a  time  delay  exists  between 
the  fall  in  the  inventory  level  and  the  fall  in  the  production  rate.  This 
time  delay  is  quite  natural  in  any  practical  situation  and  is  very  well 
presented  in  the  model  stated  above. 

The  technique  illustrated  above  could  easily  be  extended  to  optimize 
stagewise  processes  such  as  a  transportation  problem.  It  could  also  be 
extended  to  optimize  multiproduct  multi facility  scheduling  problems  with 
complex  interconnections. 
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The  choice  of  a  proper  value  for  A<f>  is  of  great  importance  if  this 
technique  is  to  be  applied  to  actual  optimization  situations.  It  has 
been  found  from  experience  that  the  value  of  A<j>  should  be  reduced  when 
the  gradient  direction  for  the  control  changes  sign  together  with  a  major 
drop  in  the  profit  function.  If  some  effective  logic  is  developed,  to 
automatically  adjust  A$  once  the  gradient  direction  changes  sign,  the 
required  computation  time  can  be  greatly  reduced.  The  author  used  a 
logic  in  which  the  value  of  Aif  was  reduced  by  half  once  the  gradients  at 
all  the  computed  grid  points  changed  sign.  This  method  failed  because 
after  the  reduction  the  value  Of  A<(>  was  either  still  too  large,  causing 
a  further  reduction  in  its  value,  or  was  too  small,  causing  the  conver- 
gence rate  to  be  extremely  slow. 

In  the  case  of  the  production  and  inventory  model  with  fixed  end 
conditions ,  it  was  found  that  by  using  values  for  the  control  variable 
far  from  the  optimal,  the  cost  C  went  above  reasonable  limits.  An 
emperical  rule  was  adopted  of  asking  for  only  part  improvement  in  the 
inventory  at  each  iteration.  Encouraging  results  were  obtained. 

As  the  problem  of  slow  convergence  near  the  optimal  still  persists, 
it  is  suggested  that  the  first  variation  of  the  gradient  technique  pre- 
sented above  should  be  used  to  get  good  starting  values  for  the  controls 
and  other  iterative  optimization  techniques  such  as  the  second  variation 
or  the  quasilinearization  technique  be  used  to  reach  the  optimum. 
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APPENDIX  A  Numerical  Solution  of  Differential  Equations 

As  seen  in  the  previous  sections,  the  gradient  method  requires  the 
solution  of  sets  of  first  order  differential  equations.   It  would  be 
appropriate  at  this  stage  to  discuss  some  of  the  numerical  techniques 
for  the  solution  of  such  differential  equations. 

Ordinary  differential  equations  are  generally  classified  according 
to  the  degree  of  difficulty  in  obtaining  numeric  solutions.  They  could 
be  classified  as 

1.  Initial  value  problems 

a)  Linear  differential  equations 

b)  Nonlinear  differential  equations 

2.  Boundary  value  problems 

a)  Linear  differential  equations 

b)  Nonlinear  differential  equations 

As  long  as  the  differential  equations  are  initial  value  problems, 
numerical  techniques  can  solve  both  linear  and  nonlinear  problems  with 
equal  ease.  However,  for  boundary  value  problems,  the  degree  of 
difficulty  in  solving  a  nonlinear  equations  is  far  greater  than  that  for 
solving  a  linear  equation. 

The  three  main  methods  generally  used  for  the  numerical  solution 
of  differential  equations  can  be  classified  as 

1.  The  single  step  technique 

2.  The  multiple  step  technique 

Under  the  single  step  technique,  the  Euler  and  Bange-Kutta  integration 
methods  would  be  discussed  while  the  Milnes  method  would  be  discussed  in 
the  multiple  step  technique. 


00 


Euler  Method  [8] 

Suppose  that  one  wishes  to  solve  the  set  of  ordinary  first  order 

differential  equations 

dx. 

4^-6^,  x2,  ....  xn,  rv   y2,  ....  ys);   i  -  1,  2,  ....  n 

J  =  1,  2 s 

where  x  is  the  state  variable  and  y  is  the  control  variable. 

The  Euler  method  is  basically  an  approximation  of  the  above  dif- 
ferential equation  in  the  form 

Ax.    dx. 

1t-"dt  =  gi  (V  V  ••••  V  yr  y2>  ••••  ys>  • 

This  approximation  leads  rather  naturally  to  the  equation 

k+1    k  _  ..     ,kk       kkk       k, 

x.    =  x.  +  At  •  g. (x, ,  x„,  ....  x  ,  y, ,  y„,  . ...  y  ) 

1      l        °i  1'   2'     '  n'  *1*  *a        s 

where  x.   indicates  the  value  of  x.  at  (k+l)At,  x.  is  the  value  of  x.  at 
i  11  l 

k   k       kkk       k 
kAt  and  g.(x  ,  x„,  ...,  x  ,  y  ,  y  ,  ...,  y  )  indicates  the  value  of  g. 

evaluated  at  time  kAt.  One  can  easily  see  how  the  above  equation  can  be 

used  in  conjunction  with  the  values  of  x.  at  zero  time  to  obtain  values 

of  x.  for  any  specified  k. 

Since  the  approximation  on  which  the  Euler  method  is  based  becomes 

exact  only  when  At  +  0  and  finite  At's  are  required  for  numerical 

calculations,  one  cannot,  in  general,  expect  the  Euler  method  to  be  very 

accurate.  Various  modifications  of  the  basic  Euler  method  are  available 

which  reduce  the  accuracy  problem. 


a? 


In  any  computation  involving  the  Euler  method,  one  must  select  the 
time  interval  At  very  carefully.  Large  At  values  lead  to  gross  errors 
vhile  small  At  values  cause  excessively  long  computation  times.  The 
suitability  of  the  Euler  method  for  use  with  a  specific  problem  depends 
upon  there  being  a  At  which  is  an  adequate  compromise  between  the  two 
effects. 

Runge-Kutta  Method 

Like  the  Euler  method,  the  Runge-Kutta  method  is  designed  for  use 
with  sets  of  first  order  ordinary  differential  equations. 

The  basic  method  is  to  write  the  Taylor  series  expansion  for  snail 
purturbation  of  the  variables  about  the  initial  conditions.  The  Taylor 
series  is  terminated  after  a  suitable  number  of  terms  and  a  series  of 
algebraic  and  operator  manupulations  is  performed  which  leads  to  lumping 
the  various  derivatives  into  terms  which  may  be  evaluated  from  formulas. 
The  most  popular  version  of  the  Runge-Kutta  is  the  system  which  results 
from  retention  of  terms  of  up  to  fourth  order  in  the  original  Taylor 
series. 

The  solution  for  a  system  of  n  equations  of  the  form 


dx. 

-£=  gi(x1,  x2,  ....  xj;   i  =  1,  2,  ....  n 


can  be  lead  through  application  of  the  formulas 


k+1    k 
x.    =  x.  +  (R.  ,  +  2  R.  „  +  2R.  ,  +  R.  , ) 
1      1     1,1      1,2     1,3    l,U 


Bt,l  =  "  •  gi(V  X2'  ••••  \] 


■i.a"4t,8i(4*fBi,l,^  +  tB2.i»  •■■\+2\,i 


So 


Si,3  =  At    •    giUl+2Rl,2;   X2  +  ?R2,2; 


k  .  1 


+  i\,2] 


and 


R. >u  =  it  •  gi(x1  +  R1)3;  x2  ♦  R2)3 


xK  +  R  J 
n    n,3 


The  index  i  goes  from  1  to  n.  As  usual,  x.  indicates  the  value  of  x^  at 
a  time  of  kit.  R.   indicates  the  J   Runge-Kutta  coefficient  for 
equation  i. 

The  use  of  these  equations  to  numerically  solve  the  system  of 
equations  is  straight-forward.  The  procedure  is: 

k   k       k 
1.  From  a  knowledge  of  x  ,  x  ,  ...,  xn>  calculate  R1  ^  R2>1; 


....  R  .. 

n,l 

k   k 
2.   From  x  ,  x  , 

Rl,2;  R2,2; 


,,  x  and  R.,  ,;  R„  .;...,  R,  calculate 
n      1,1   2,1       n,l 


n,2 


3.  From  xx,  x£,  . . . ,  xr  and  R   ;  ^2,2'' 

R 


Rl,3;  R2,3; 


k   k 


n,3' 


h.     From  xf,  XT,  ....  X  and  R   ;  R 


1'  2 
R1,U;  R2,l*; 


1,3*  "2,3' 


. ,  R  .,  calculate 
n,2 


i  ,  R  „ ,  calculate 
n,3 


,,  R 


n,V 


5.  Calculate  x.   from  x^  Ri  ^  RA  2;  R^^  3  and  R±  j^. 

6.  Repeat  steps  1  ■*  5  until  the  final  value  of  A  is  reached. 


The  problem  of  accuracy  is  not  of  great  importance  in  the  Runge-Kutta 
integration  scheme  as  the  truncation  error  is  of  the  fifth  order.  However, 
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the  stability  problem  sometimes  arises  as  the  process  might  not  converge. 

Multiple  Step  Method  [8] 

For  the  solution  of  the  differential  equation  of  the  type 

£.*(..*>. 

the  formulas  for  the  recurrence  relationship  could  be  represented  by 

*(W  -  x(tk-r' 

=  h(  t,  f(x(tk),  tk), 

f(x(Vl),  xCt^),  ...  f(x(tk_n),  tk_n)) 

where  r  and  n  are  positive  integers. 

To  evaluate  x(tk+1),  the  values  of  x(tk),  *(*fc_i)>  ■••«  x^k-n' 

and  x(t   )  must  be  known.  Hence  it  is  seen  that  it  is  not  possible  to 
k-r 

calculate  x(t   )  directly  from  the  initial  value  x  .  Also,  to  start 
the  calculation  the  points  x(t  ,),  x'tk-2^'  -"  must  be  kn0WI1  throuSh 
another  integration  method. 

To  increase  the  accuracy,  two  integration  formulas  are  generally  used 
in  the  multiple  step  method.  The  first  formula,  known  as  the  open  end 
integration  formula,  is  used  to  predict  the  approximate  value  of  xt^.^)' 
Then  the  second  formula,  or  closed  end  formula,  is  used  to  generate 
a  more  accurate  x(t,  .).  This  latter  formula  may  be  iterated  to  obtain 
as  accurate  an  answer  as  desired. 

These  two  formulas  form  a  predictor  correcter  scheme  which  is  a 
powerful  numerical  tool.  Milnes  method  is  probably  the  best  known 
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multiple  step  integration  formula.  The  predictor  for  this  method  is 


x(W  ■  x(\-3)  +  3  At  [2f(x(V-  V  -  fWVi>>  Vr 


+  2f(x(tk_2),  V2))] 


and  the  correcter  is 


To  begin  the  integration,  the  starting  values  at  the  three  grid 
points  t  ,  t    and  t    can  be  obtained  by  a  single  step  integration 
formula  or  by  using  Taylor  series. 

The  single  step  methods  have  a  number  of  advantages  in  terms  of  the 
use  of  digital  computers.  First  in  using  the  multiple  step  methods  the 
starting  values  must  be  calculated  by  some  other  methods;  no  such 
predictions  or  corrections  are  necessary  for  the  single  step  methods. 
Second,  during  the  integration  process  several  different  values  of  the 
integration  step  At  may  be  necessary  in  solving  the  sane  equations.  It 
is  not  easy  to  reduce  the  integration  step  At  for  the  multiple  step 
methods  as  the  integration  proceeds.  Some  kind  of  interpolation  formula 
must  be  used  to  reduce  this  step  size. 

Because  of  its  high  relative  accuracy  and  ease  of  computation,  the 
Runge-Kutta  method  is  used  for  the  numerical  solution  of  the  differential 
equations  encountered  in  this  report. 
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ATPKNDIX  S3  Plow  Ghnrta 


Start 


Dimension  pjlcl 
Function  Statements 


Calculate  partial 
derivatives  of  final 
perfornEiice   criterion 
with  respect   to   control 
variables 


H 


Calculate  a  nev   control 
sequence  usin£  old  control 
sequence,  partial  derivetives 
and  step  size  Delta.   


Head  initial  conditions 
end  parameters 


Asgu.'G  a  control 
seauence 


Calculate  state  voxiables 
and  final  value  of  objective 
function  corrospondirv;  to 
assurred  control  sequence 


Calculate  state  variables 
and  objective  function 
corroopondinc  to  new  control 
variables. 
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APPENDIX  C 
COMPUTER  PROGRAKS 


)    F0RMAU8F8.3,  11,13)  qP- 

.  F0RMATIF8.6)  ^ 

)    F0RMATI1H  ,  10HITERATIOM  ,13) 

)  FORMATdH  ,3f.l6.8) 

)  FORMAT! 1H  ,51H  SUM    SUM  OF  DTHETA/DZSUM  OF  DTH/DZ+DS/OZ I 

DIMENSinNT(200),Pl(200),P2(200) , P3I200) , P4I2C0) .QK2CO) ,02(200) , 
103(200) ,0^(200 ), St  200) f XI ( 200 ) , DX ( 200 ) , DX2 ( 200 ) , X2 ( 200) .SXK200) , 
2ST1200) , E 1200) ,  Z (200) ,PSZ( 200) 

DO  10  1=1,101 

H=l-1 

X=H*D 
I  T( I )=X*4.75+3. 

I  SUM=0 
SUMT=0 
SUMP=0 

DO  A  1=1, 101 

S  ( I  1  =  1-1 

Xl(l)=5. 

X2(l)=0. 

P1(I)=(T( t)-A-B«D«S(I) )*D 

P2( I )  =  D*( T( I )-A-B*(D*S( I  1+0/ 2.  )  ) 

P3(I)=P2( I) 

P4(  I  )=D*(T(I  )-A-F-*(D*S(  I  )+D)  ) 

DXI I )  =  I 1./6. )* (PI  I  I  )+2.*P2( I K2.*P3(  I )+PM I  )  ) 

X1(I+1)=X1( Il+CXI I ) 

QUI  >=(CI*1  (Xl( I)**2)+CP*EXP(  (PA-TI  [) )**2) )*D 

02(1 )=(CI*(K1(I)*.5*PI (I) )**2+CP*EXP( I PA-TI I ) )**2))*D 

Q3<I)  =  (CI*(XKI)  +  .5*P2U))**2  +  CP*EXP((PA-T(I))**2))*D 

QMI)=<CI*(X1(I)+P3U) )**2+CP*EXP( ( PA- T ( I ) ) **2 ) ) *D 

DX2II  )  =  ( 1./6. )*(QltI )+2.*Q2( I)+2.*Q3 

DX2(I)  =  (l./6. )  MQK I  )+2.*Q2( I)+2.*Q3(I)+Q4(  I)) 

X2( I  +  1)  =  X2( I  1+0X21 t  I 

E(I)=S( [)*0 

CONTINUE 

FORMAT! 1H  ,.« I T£R« *3X, »TIME» ,2X, ' INVENTORY' , 4X,' COST ',  10X, ■ THE TA' , 

II  OX, 'CS/DX1',1CX,'DS/DZ',10X, 'DTH/DZ' ) 
PRINT  300 

DO  5  1=1,100 

SX1(1CD  =  0. 

N= 101-1 

SXl(N)=SXl(Ntl)+2.*D*CI*XllN) 

ST!N)=0*(SX1(N+1)-2.*CP*(PA-T(N) )*EXP( (PA-T(N) )**2)  ) 

SUM=SUM+ST(N)**2 

ST(101)=ST!100) 

SUM=SUM+ST( 1011**2 

DO  15  1=1,101 

TZ(I)=D 

sumt=suht*:tz(  i  1**2 

00  6  1=1,101 

PRINT  200,  I,  El  I  ),X1(  I),X2(  I),T(I),SXK  I  ),ST(I),TZ!I  ) 

PRINT  700 

PRINT  500, SUM, SUMT, SUMP 

K=K+1 

IF(K-L)7O,50,50 

0TH=9.25-X1(101) 

A2=(DEL*SUMP-DTH*SUM)/(SUMP**2-SUMT*SUM) 

A1=(DEL-A2*SUMP)/SUM 

DO  7  1=1,101 

T(I)  =  T(I)<-(A1*ST(I)  +  A2*TZ(  I  )  ) 

PRINT  ' :0,K 


DIKENSI 
DIMENSI 
FORMAT! 

.2X.13HA 
FORMAT! 
FORMAT! 
FORMAT! 
FORMAT! 
FORMAT! 
FORMAT! 
FORMAT! 
FORMAT! 
FORMAT! 
FORMAT! 
READ  18 
READ  14 
READ15, 
PRINT19 
PRINT16 
PRINT17 
DO  9  1  = 
AII)=2. 
SUM=0 
DO  1  1  = 
T=I-1 
T=T*D 
P(I)  =  AB 
A1=(P(I 
B1=!(C+ 
A2=(P!I 
B2=(!C+ 
B3=(  (C+i 
B4=(  (C+ 
A3=(P(I 
A4=(P(I 
C1=(F*X 
C2=(F*( 

L)*(X2(I 
C3=(F*( 

M*(X2(I 
C4=(F*( 

L1+B3) )* 
X1II+1) 
X2II  +  1) 
X3I1+1) 
E(I)=T 
DO  2  1  = 
SXK101 
SX2I101 
5X3(101 
N=101-I 
SX1(N)= 
SX2(N>= 

L-SX3IN  + 
SX3(N)= 
ST(N)=( 

L*D 

SUM=SUM 
ST(lOl) 
SUM=SUN 


OMP!  I 
0MSX1 
1H  ,"> 
DVE'U 
IH  ,F 
1H  ,9 
E 1  5  .  9 
8F8.3 
2F5.2 
IH  ,8 
IH  ,2 
3F10. 
1 H  ,3 
F20.B 
,  X 1 1  1 
,AB,R 
D.DEL 
,XL  !  1 
,AB,B 
iD.DE 
1,101 
5 

1,101 


02)  ,X 

!  1C2) 

II    TIM 

I  SEME 

5  .  2  ,  ft 

HITFR 

) 

,214) 

) 

F8.3f 

F5.2) 

2) 

F10.2 

) 

),X2( 

,C,PO 

)  ,X2( 
,C,PO 

L 


2!  102) , XI (102) ,  X3I102) ,A( 102) ,E( 102) 

,SX2(120) ,SX3(120) ,ST(102) 

E,6X,9H INVENTORY, 10X , 5HSALES , 9X, 6HPR0F I T , 

NT,5X,10HPR0DUCTI0N,7X,8HGRADIENT) 

C15.4) 

AT  ION, 16) 
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214) 


1),X3!1) 

,PI  ,F,CI ,CA,K,L 


1) ,X3(1) 

,PI ,F,CI ,CA,K,L 


+  B*T 

)-X2( 

A(l)  ) 

)-X2< 

AID) 

All.)  ) 

4(1)1 

1-X2I 

)-X2( 

2!  I)- 

X?(I) 

1+B1/ 

X2(I) 

)+B2/ 

X2II) 

D 

=  X1(  I 

=  X2<  I 

«X3*(I 

1,100 

1  =  0 
)=0 

)  =  1 


1)1*0 

*(X2(l)-IX2( I)**2)/P0) )*D 

I  1-B1/2. )*D 

*(X2(I)+Bl/2.-((X2(IH-Bl/2.)**2)/P0))*D 

*(X2(  t  HB2/2.-I  1X2  1  I 1+B2/2. )**2)/P0) 1*0 

*(X2[ I 1+B3-! (X2( I ) +R3 ) **2 ) / PO ) )*D 

()-B2/2.)*D 

I  l-B3J*0 

( (PI-X1  (  I  ) )**2)*CI-CA*(A( I )**2)*X2( I)  )*D 

+B1/2.1-I (PI-IXK I 1+A1/2. ) )**2)*CI-CA*!A! I  1**2 

2.)  J*0 

+B2/2. )-( (PI-IXK I 1+A2/2.) ) **2 ) *C I-CA* ( A ( 11**2 

2.  )  1*0 

+B3)-((PI-(X1(I)+A3) )**2)*CI-CA*IA! I )**2)*(X2( I 

l  +  (AH-2.*A2  +  2.*A3  +  A4)*l./6. 
)+<Bl+2.*B2+2.*B3+B4)*l./6. 
)+(Cl+2.*C2+2.*C3+C4)*l./6. 


(SX3!N*l)*2.*CI*!PI-Xl(N)))*D 

l-SXKN+1 )+SX2(N+l)*(C+A(N) )*( I  .-2. *X2 (N ) /PO) 

IN)**2)-F) )*0 


SXKN+D  +  C 
SX2IN+1 1+ I 
1)*<CA*(AP 
SX3IN+1I 

SX2IN+1 1*1X2 (N)-(X2(N)**2)/P0)-SX3(N+1)*2.*A(N)*X2(N)*CA) 

+ST(N)**2 
=ST(100) 
♦ST! 1011**2 


PRINT10 

DO  3  1=1,101,5 

PRINT  11.  ,E([),X1([),X2(I),X3(I),A(I),P!I),ST(I) 

K=K+1 

IFJK.GE.18)DEL=.05 

IF(K-LK,5, 5 

DO  6  1=1,101 

A(I)=A( I)40£L*ST( I 1/SUM 

PRINT12.K 

GO  TO  7 

DO  8  1=1,101 

WRITEIZ, 13)M  I) 

STOP 

END 


100. 


2. 


10. 


,  IS 


1.5 


120 


99 


DIMCNSICNTK  10 
DIMENS  ION  Xld 
DIMENSICNSX2I1 
DIMENSI0NST2I1 
DIMENSI0NZX5!  1 
DIMENSIONPSZH 
DIMENSION    AA1 I 

1DB1TK  10?)  ,l)A2 
DIMENSION  CLT1 
FORMATdH  ,6HD 
FORMAT! 1H  ,2F8 
FORMAT! 1H  ,7F8 
F0RMAT11H  ,3F1 
F0RMATI1H  ,7F1 
FORMAT! 1H  ,2G2 
FORMATdH  ,2F1 
F0RMAT(3f 11.4) 
FORMAT! 1H  ,3F6 
F0RMATI1H  ,F4. 
FORMATdH  ,F4. 
FORMAT! 1H  ,.4HD 
FORMATdH    r4HT 

1C.    2A,  7X.RH 

FORMATdH    ,4HT 

1R0FIT,10X,5HGR 
FORMATdH  ,9F1 
FORMAT! IH    ,3X, 

15HSUMZ3,7X,5HP 
FORMATdH    ,9HI 
F0RMATI3E20.-3) 
F0RMAK2F4.2) 
REAC7C0,X0,Y0 
F0RMATI7F8.2.3 
READ6C0,Xl( I) , 
FCRM,AT(8F10.4) 
F0RMATI7F10.4) 
READ8C0,C0,C,P 
READ850,CI,CA, 
FCRMATI2E10.3, 
READ9C0,GA,Gn, 
F0RMATI5F6. i) 
FCRMAT12F10.2) 
READ    150,D,-DEL 
DO    11     1=1,101 
T1(I>=340 
T2( I)=340 
T3(I)=2. 
OLTK  I)  =  TK  I  ) 
0LT2! I)=T2!  I) 
0LT3! I )  =  T3(  [  ) 
PRINT10l,X0,Y0 
PRINT1Q2.X1!  1  ) 
PRINT103.CQ.C, 
PRINT104,CI,CA 
PRINTI05.GA.GB 
PRINT106.D.0EL 
DO    15    1=1,101 
PRINT108,TUI  ) 
A=10. 
SUM1=0 


2), 
C2) 
C2) 
C2) 

C2) 
102 
102 
1 2  I 
(10 
ELT 
.2) 
.?, 
0.4 
C.4 
C.  I 
C.3 


T2I102) ,TS (102), E( 102) ,S (102), SX 1(102 ),SY 1(102) 

,Yl(102),X2(102),Y2!102),X3(102),X4d02),X5d02) 

, SY2(102),SX3(102),SX4(102),SX5( 102) ,ST1< 102) 

,  ST  3  (102),  ZXK  102)  ,  ZX2(  102)  ,ZX3(  102),.  ZX4  1102) 

,ZY1(102),ZY2(102),ZT1(102),ZT2(102),ZT3(102) 

),r>SZ2(102),PSZ3(  102) 

), AO 11102) ,AA2(102),AB2! 102 ) , DA  IT  1  (  102  ) , 

102) ,DB2T2(102) 

2)  ,.GLT2(  1021.0LT3!  102) 

l=,F9.6,6HDELT2=,F9.6,6HDELT3=,F9.6) 

215) 

) 

) 

,2F10.0) 

) 


.2) 

2.6CI 

2,  3E1 

EL=,F 

IME.7 

CCNC 

IMF, 6 

,T1,1 

2.6) 

4HSUM 

SUM1, 

TFRAT 


[2] 


5.4) 

5.4.E18.7.3E15.4) 

5.2) 

X.8HC0NC.     1A,7X,8HC0NC1    B,9X,6HTEMP    1,7X,8HC0N 

2    B.9X.6HTEMP    2) 

X,9HIMVENTORY,10X,5HSALES,7X,8HADV.COST,12X,6HP 

0X,5HGR.T2,10X,5HGR. T3) 

l,8X,4HSUM2,.8X,4HSUM3,7X,5HSUMZl,7X,5HSUMZ2,7X, 

7X,5HPSUM2,7X,5HPSUM3) 

ION, 14) 


Villi  ,X2d)  ,Y2d),X3(  1),X4(1),X5(.1),K,L„M 


,AM,T1M,Z1,Z2,  Z3 

R,Q,Vl,V2,CT 

2F7.0) 

EA.EB 


Tl,DELT2,DELT3,OLCOS 


,Y1(1),X2(1) ,Y2(1),X3(1),X4(1) ,X5(1),K,L 

P,AM,T1M,Z1,Z2,Z3 

,R,Q,V1,V2,CT 

,EA,EB 
Tl 

,T2d  )  ,T3(  I  ) 


100 


SUM2=0 
SUM3=0 

SUP;Z1  =  0 

SUMZ2=0 

SUMZ3=0 

PSUM1=0 

P5UM2=0 

PSUM3=0 

SUMT1=0 

SUMT2=0 

SUyT3=0 

0011=1,101 

S  ( I )  =  1-1 

AAK  I  J=GA*E 

ABKI  )=GB*E 

AA2(I )=GA*E 

AB2( I )«GB*E 

A1=(Q*(X0-X 

A2=(0*(X0-X 

A3=(Q*(X0-X 

A4=(0*(X0-X 

Bl=(Q*lYG-Y 

B2=(G*( YO-Y 

1/2. ) )*C 
B3=(C*(Y0-Y 

1A2/2. ) )*D 
B4=(Q*(Y0-Y 
Xl(I+l)=Xl( 
Y1(I+1)=Y1( 
C1=(Q*(X1( I 
C2=(Q*(X1( I 
C3=(C*(X1( 1 
C4-(Q*U1(  I 
D1=(Q*(Y1I t 
D2=(Q*(Y1 ( I 

1(X2( I 1+C1/2 
D3=(C*(Y1( [ 

KX2IT1+C2/2 
D4=(0*(Y1<  I 

1)  )*0 
X2( l+l)=X2( 
Y2(I+1)=Y2( 
Fl=( (C+T3II 
F2=( (C+T3I I 
F3=( (C+T3I  [ 
F4*UC+:T3M 
XA(  I  +  I)  =  XM 
E1=(C*Y2(I) 
E2=(Q*(Y2( I 
E3=(Q*(Y2( I 
E4=(Q*(Y2(  I. 
X3(I+1)=X3( 
G1=(Z1*CQ*X 

1-CT*(  (TIM-T 
G2=(Z1*CC*( 

1(  11-01/2. )- 

2ITHI  )-T2(I 
G3=(Zl*CQ*l 

KT1-D2/2.  )- 

2(T1( t)-r2(  I 


XP<-EA/(R*T1(I  )  )  ) 

XP(-EB/(R*T1 ( I  )  )  ) 

XP(-EA/(R*T2( I )  )  ) 

XP(-E3/(R*T2( I ) ) ) 

III) l/Vl-AAK I )*X1( I ) )*D 

l(I)-Al/2.)/Vl-AAKI)*(Xl(I)+Al/2.))*D 

l( I  1-A2/2. 1/V1-AA1 ( [)*(X1( D+A2/2. ) ) *D 

11  I )-A3) /Vl-AAK I)*(X1(I)+A3))*D 

1(1)  )/Vl-AEl(  I  )*Y1(  IHAAK  I  >*X1I  I)  )*D 

l(I)-Bl/2.)/Vl-ABl(I)*(Yl(I)+31/2.)+A41(I)*(XKI)+Al 

1(1  1-D2/2. 1/V1-ABK I)*(Y1( I)+B2/2.)+AAU 1 )*(Xlt  I  )  + 


llll-B 
Iltll. 

n  +  ii. 

)-X2( I 
I-X21I 
)-X2(  I 
I-X2I  I 
)-Y2 ( I 
)-Y2(I 
.1  )*C 
)-Y2t I 
.» J*C 
)-Y2(I 


31 /Vl-ABK I)*(Y1(I)+B3)+AA1(I)*(XK I1+A3) )*0 
/6.)*IAl+2.*A2+2.*A3+A4) 

/6.)*(Bl+:2.*B2  +  2.*B3  +  B4) 
) J/V2-AA2I I )*X2I1) )*C 

)+Al/2.-Cl/2.)/V2-AA2(I)*(X2(I)+Cl/2.) 1*0 
1+A2/2.-C2/2. 1/V2-AA2I I ) *  I X2 ( I > +  C2/2 . )  ) *D 
)+A3-C3)/V2-AA2( I ) * ( X2 ( I ) +C3 ) ) *D 
))/V2-AB2(I)*Y2(I)+AA2(I)*X2U))*D 
H-B1/2.-D1/2. 1/V2-AB2I I)*(Y2(I)+01/2.)+AA2(I)* 

1+B2/2.-D2/2. I/V2-AB2I I ) * ( Y2 ( I ) +C2/2 . ) +A A2 ( I ) * 

)+B3-DJ)/V2-AB2( I ) * ( Y2 ( I ) +03 ) +AA2 ( I ) * ( X2 ( I ) +C3 


1  !  t 
I  )  ► 
)  I* 
!  M 
)  »* 
))* 
I  )  + 
-CC 
)+D 
I  t-0 
]  fD 

I)  + 

M  I 

I  ( I 

XM 
CI* 
)  )* 
X4( 
CI* 
)  )* 


(1./ 

(1./ 

XM  I 

(X4( 

(X4( 

(X4( 

(  I./ 

*X4( 

1/2. 

2/2. 

3)-C 

(1./ 

)  +  Z2 

))** 

I  )+F 

(AM 

*2) 

I  >  +  F 

(AM 

♦  2)- 


6.  ) 
6.  ! 
)*( 
I  )  + 
I  )  + 
I  >  + 
ft.  ) 
I  )  ) 
l-C 
l-C 
Q*( 
6.  ) 
*Q* 
2+1 
1/2 
X3( 
CA* 
2/2 
X3( 
CA* 


*(C1+2.*C2+2.*C3+C4) 

*(D1+2.*D2+2.*D3+D4) 

1.-X4I Il/P) 1*0 

Fl/2. )*U.-IX4( I l+Fl/2. )/P))*D 

F2/2. )*ll.-(X4( I)+F2/2.)/P))*Q 

F3)*(l.-(X4l I 1+F31/P) )*D 

*(Fl+2.*F2+2.*F3+F4> 

*D 

G*(XM  I  )+:Fl/2.  )  )*D 

G*(X4( I I+F2/2.) )*0 

X4( I )+F3) )*D 

*(F1+2.*E2+2.*E3+E4) 

X2( [ )+Z3*G*(l.-X2( I (-Y2II ) )-C I* ( AM-X3 ( I 

Tl( l)-T2(I) )**2)-CA*( (r3(I)*X4(I))**2.) 

,)  +  Z2*Q*(X2(  D+Cl/2.  >+Z3*Q*(  1.-X2I  I  )-Cl 

I)-El/2.)**2  -CT*( IT1M-T1II1 

(  (T3( I)*(X4(I)+Fl/2.))**2.>  >*D 

. )+Z2*Q*(X2( D+C2/2. )+Z3*Q*( 1.-X2I I )-C2 

I)-E2/2.)**2  -CT*(  (T1M-TKI) 

(  (T3(I  )*(XM  D+F2/2.)  )**2.)  )  *0 


)  )**2 
)*D 

/2.-Y2 
1**2  + 

/2.-Y2 
)**2  + 


101 


G4=(Z 
1-CI*( 
2*2)-C 

X511  + 

E<n* 

DO    2 

SXK1 

SYK1 

SX2I1 

SY2(  1 

SX3I1 

SXM1 

SX5I1 

N=101 

SXKN 

1*Q/V2 
SYKN 
SX2IN 

11)*(Z 
SY2IN 

13*Q)) 
SX3IN 
SX4IN 

1-SX5I 
SX5IN 
DA1T1 
DB1T1 
DA2T2 
DE2T2 
STUN 

1(N)*D 
ST2IN 

1*DB2T 
ST3(N 

IN) )*D 
SUMT1 
SUMT2 
SUMT3 
SUM1  = 
SUM2  = 
SUM3  = 
STK1 
ST2(  1 
ST3(1 
SUM1  = 
SUM2  = 
SUM3= 
SUMT1 
SUMT2 
SUMT3 
IFIX5 
00  37 
Till) 
T2(I) 
T3(I) 
IF  (  SO 
IFIOL 
GO  TO 
IFIOL 
IFISU 


1*CQ*(X4 
AM-X3I [ ) 
A*( (T3( I 
l)*X5U) 
S(I)*D 
1=1 , LOO 
0 1 )  =  0 

on=o 

01)=0 

ou=o 

011  =  0 
011*0 
01)  =  L 
-I 

i  =  SXUN  + 
)*D 

)=SYL(M+ 
)=SX2(N+ 
2*Q-Z3*0 
)=SY2(M+ 
*C 

)=SX3(N+ 
)  =  SXMM+ 
N+l)*(2 
)  =  S  X  5  (  I  + 
IN)*(5A* 
(N)«IGB* 
<N)=(GA* 
1N)=(GQ* 
)«(SX1(N 
B1TUM)  ) 
)=(-SX2( 
2(N) ) fSX 
)=(SX't(M 

■sumti+:s 

=SUKT2*S 
=  SUMT3+'S 

sumi+:sti 

suM2+:sr2 

SUM3+ST3 
01)=ST1( 
01)=ST2< 

cu  =  sr3( 
suMi+:sri 

SUM2+ST2 

suM3+:sr3 

=SUMT1+S 
=SUMT2+S 
■SUHT3+;S 
( 101). Gc 

1=1,101 
=0LTL( I  1 
=  0LT2(  I) 
=CLT3( [  ) 
MT1.GT.0 
ST1.GT.0 

31 
ST1.LT.0 
NT2.GT.0 


(I)+F3)+Z2*0*(X2( I  )+C3) 

-E3)**2  -CT* 

) *{X4( I )+F3) 1**2. ) )*D 

+  (  l./6.)*(Gl+2.*G2+2.*G3+G41 


+Z3*Q*(  1.-X2I D-C3-Y2I I1-D31 
( (TlH-Tl(I) 1**2+ (Tit I1-T2I I ))* 


1)-MSX1M+1)*(-Q/V1-<\A1(N)  )  +  SY1  (  N+  1  )  *AA  1  (N ) +SX2 ( N+ 1 1 


.IN) )+SY2(N+l)*Q/V2)  *D 

!(N) )+SY2(M+l)*AA2(N)+SX5IN+ 


1)+(SY1(N+1)*(-C/V1-AB1I 

1  )  +  (SX2(N+l)*(-Q/V2-AA2( 

)  )*D 

1)+(SY2(N+1)*(-G/V2-AB2(N) )+SX3(N+l)*Q+SX5(N+l)*(-Z 


1 1+SX5 (N  +  l )*2.*CI*(AM-X 

l)  +  (-SX3(sJ  +  l)*CQ+'SX4(N  + 

*CA*{T3(N)**2l*X4{N) 

1) 

EA/(R*(T1(N»**2) ) )*EXP( 

EB/U*(T1(N)**2) ) )*EX?( 

E«/(R*( T2(N)**2) ) )*EXP( 

EB/(R*(T2(N)**2) ) )*EXPt 

+1)*(-X1 (N)*0A1T1 (N) )+S 

+SX5(N+1)*(-CT*(-2.*(T1 

NU)*X2(N)*D42T2(N)  +  SY2 

5(\+l)*(CT*(2.*(Tl(N)-T 

*l)*X4tN)*U.-X4(N)/P}- 

T 1  (  N ) 

T  2  (  M  ) 

T3(N) 

(M**2 

(M**2 

(\')**2 

ICO) 

ICC) 

ICC) 

1101 )**2 

( 1011**2 

( 1011**2 

Tl( 101 ) 

T2I101  ) 

T3(L01) 

•CLCOSIGO  TO  113 


)G0  TO  30 
. )CELTl=0ELTl/2. 


3(N) )*D 

1)*(C+T3(N) )*(1.-2.*X4(N)/P) 

-CQ*Z1) )*D 

-EA/(R*T1(M) ) ) 
-EB/(R*T1(N) ) ) 

-EA/(R*T2(N) ) ) 
-E3/(R*T2(N) I ) 
YKN+1)*(+X1  (N)*CA1T1(N)-Y1 
M-TKN)  )+2.*(Tl(N)-T2(N)  )  )  )  )*C 
(N+l)*(X2(N)*DA2T2(N)-Y2(N) 
2(N))J))*0 
SX5(N+l)*2.*CA*(X4(Nl**2)*T3( 


)CFLTl=DELTl/2. 
)G0  TO  32 


102 


IFIOLST 
GO  TO  3 
IFIOLST 
IFISUMT 
IFIOLST 
GO  TO  3 
IFIOLST 
PRINT36 
GO  TO  2 
0LST1=S 
0LST2=S 
0LST3=S 
DO  3  1  = 
ZXK101 
ZYK101 
ZX2I101 
ZY2( 101 
ZX3I101 
ZX4I  101 
ZX5I101 
N=101-I 
ZX1(N)= 

L*Q/V2)* 
ZY1(N)= 
ZX2(N)= 

tl)*(Z2* 
ZY2(N)= 

L3*Q) )*0 
ZX3(N)= 
ZX4(N)= 

1-ZX5IN  + 
ZX5<N)= 
ZTKN)* 

KN)*DB1 
ZT2(N)= 

L*DB2T2( 
ZT3(N)= 

LN))*D 
SUMZ1  =  S 
SUMZ2=S 
SUMZ3=S 
ZTK101 
ZT2I101 
ZT3I101 
SUMZ1=S 
SUMZ2=S 
SUMZ3=S 
DO  4  1  = 
PSZ1I  I) 
PSZ2I  I) 
PSZ3I  I) 
PSUM1=P 
PSUK2=P 
PSUM3=P 
IF(M)21 
PRINT1C 
DO  5  I 
PRINT20 
PRINT  1 
DO    51     I 


2.GT.0 

3 

2  .  L  T .  0 

3.GIV0 

3.GT.0 

5 

3.LT.0 

.DELTl 

0 

UMT1 

UKT2 

UMT3 

1,1.00 

)  =  0 

)=0 

1=0 

)  =  0 

)  =  1 

)  =  0 
)*0 


)CELT2=DELT2/2. 

.)CELT2=DELT2/2. 

. IGO    TO    34 

. )CELT3=DELT3/2. 

. >0ELT3=DELT3/2. 
,CELT2,DELT3 


ZXlH-tl  )  +  (ZXl(N  +  l)*(-Q/Vl-AAl(N)  )  +  ZY  I  (  N+  1 )  *AA1  (  N  )  +  ZX2  (  N+l ) 

D 

ZYLI  J+l 

Z  X  2  I  M  +  ] 

Q-Z3*Q) 

ZY2IN+1 


1)+(ZY1(N+1)*I-Q/V1-AB1(N) ) +ZY2 I N+ 1) *Q/V2 >     *D 
1)+IZX2(N+1)*(-Q/V2-AA2(N) )<-ZY2(N+l)*AA2(N)+ZX5(N  + 

!)  )*D 

tlt  +  (ZY2(Ntl)*(-Q/V2-AB2(N)  ) +ZX  3  I  N+ 1 )  *Q  +  ZX5  I  \l  + 1 )  *  (-Z 


ZX3(i>l*-l)+ZX5(N+l)*2.*CI*(AM-X3(M)  )  *0 

ZX4{N+1 )♦  I-ZX3  1N  +  1 )*CQ+ZX4(N+1)*(C+T3(N1 )*( l.-2.*X4(N)/P) 

l)*(2.*CA*(T3(N)**2)*X4(N)  -CQ*Z1) )*D 

ZX5IN+1) 

IZX1(N+L)*(-X1CN)*DAIT1(N) ) +ZY1 ( M+ t ) * ( +X1 ( N ) *DA 1T1 (N )-Yl 

TKN)  )+ZX5(M+l)*(-CT*(-2.*(TlM-Tl(N)  )+2.*(THN)-T2(N)))))*D 

(-ZX2I  kJ  +  l)*X2(N)*DA2r2(N)  +  ZY2(N+l)*(X2(N)*CA2T2(N)-Y2(N) 

M)  )+'ZX5(N  +  l)*(CT*(2.*ITl(M)-T2(N>  )  )  )  >*D 

(ZX4(N+1)*X4(N)*(1.-X4(N)/P)-    ZX5 ( N+ 1 ) *2 . *CA* I X4 ( N ) **2 ) *T3 ( 


UVZI+ZT1(N)**2 
UMZ2+:ZT2(N)**2 
UMZi+'ZT3IN)**2 

)  =  zrn  ico 

)=ZT2(1G0) 
)=ZT3(1CC) 
UMZH-ZU(l0l>**2 
UMZ2+ZT2U01)**2 

UMZ3+:ZT3(101  1**2 

1,101 

=ST1( I )*ZT1( I) 

=  ST2( I )*ZT2( I ) 

=ST3( I )*ZT3( I ) 

SUMH-PSZll  I  ) 

SUM2<-PSZ2(  I  ) 

SU«3+'nSZ3(  I  ) 

,21,22 

0 

1,101 

0,E(I),X1([),Y1(I),T1(I),X2(I),Y2(I),T2II) 

25 

=1,101 


io3 


PRINT  250,E(I),X3<I)  ,X4(I)  ,T3(I),X5(I),ST1II),ST2(I  ),ST3(  1) 

GC  TO  111 

PRINT  100 

DO  56  1  =  1, 101 t 5 

PRINT200,E( Il.XK I )  ,YK I)  ,T1(  I),X2(  I) ,Y2(  I) , T2 ( I) 

PRINT    125 

DO    57    1=1,101,5 

PRINT    2  50.EU ) ,X3(I ) ,X41I ) ,T3( I),X5(I >,ST1 I  I), ST 21  I  ),ST3(I) 

PRINT  3C0 

PRINT  400,SUM1,SUM2,SUM3,SUMZ1,SUMZ2,SUMZ3,PSUM1,PSUM2,PSUH3 

K=K+1 

IF(K-U6,10,10 

DELZ=10.-X3( 101) 

TSUMl  =  SUMl+:SUr2  +  SUM3 

TPSUM=PSUK1+PSIN2+PSUF3 

TSUMZ=SUMZl+SUf Z2+SUMZ3 

A1=DFL/TSUM1 

0LC0S=X5( 101) 

DO  7  1=1,101 

OLTK  I  )  =  T  1  (  I) 

0LT2( I)=T2<  I) 

0LT3(I)=r3( I) 

Tl(  I  )=T1(  I  ltDELTl*       STKD/SUM1 

T2( I)=T2< I X-DELT2*       ST2( I 1/SUM2 

T3U)=T3( I )+DELT3*ST3( I 1/SUM3 

IFIT3I  n.LT.O.  )G0    TO    71 

GO    TO    7 

T3(I 1=0 

CONTINUE 

PRINT    5C0.K 

GO    TO    20 

DO    61    1=1,101 

WRITE(2,251)     mill  T?(  I  ),T3(  I) 


STOP 

END 
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The  objective  of  this  report  is  to  apply  the  functional  gradient 
technique  to  optimize  management  systems.  The  basic  methodology  in  the 
functional  gradient  technique  is  to  obtain  a  set  of  functional  equations 
which  gives  the  gradient  of  the  objective  function  with  respect  to  the 
control.  The  objective  function  is  then  improved  in  this  gradient 
direction.  If  the  state  variables  have  to  satisfy  certain  final  conditions, 
a  penalty  function  is  introduced  in  the  equations. 

A  set  of  problems  in  the  field  of  production,  inventory  control  and 
advertising  have  been  solved  with  the  purpose  of  making  a  critical 
study  of  the  technique.  The  first  problem  solved  considers  sales  to  be 
fixed  and  production  controlled  in  order  to  minimize  cost  and  maintain 
a.   certain  desired  inventory  level  at  the  final  time.  In  Section  3.2 
the  diffusion  model  is  used  to  determine  sales  with  the  production 
considered  as  a  constant.  The  control  here  is  advertising  and  the  objective 
is  to  maximize  profit.  In  Section  3.3  a  set  of  three  problems,  each  with 
six  state  variables  and  three  control  variables,  is  solved.  The  results 
obtained  are  satisfactory  but  the  inherent  difficulty  of  slow  convergence 
rate  near  the  optimal  still  persists. 


